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Linear  Stability  Analysis  of  Hypersonic  Boundary  Layers 

by 

Srinivas  Tadepalli  j  and  Joel  H.  Ferziger  t 
Stanford  University 
Stanford,  California  94305. 


1.  Abstract 


For  linear  stability  analysis  of  high  speed  boundary  layers,  we  use  a  direct  method  for 
its  ability  to  yield  eigenvalues  without  a  priori  knowledge  and  to  capture  all  modes. 
Temporal  linear  stability  analysis  is  performed  for  the  2D  boundary  layer  on  a  flat 
plate  using  the  local  parallel  flow  assumption.  An  estimate  of  all  the  eigenvalues 
is  obtained  by  solving  the  generalized  eigenvalue  problem  (Malik,  1986);  a  local 
eigenvalue  search  is  used  to  improve  the  accuracy  of  the  most  unstable  eigenvalue. 
A  compact  fourth  order  accurate  method  is  to  compute  both  the  mean  flow  and 
the  most  unstable  eigenfunction.  This  method  is  more  efficient  than  lower  order 
methods.  Local  grid  refinement  based  on  error  estimates  is  useful  in  providing  the 
accuracy  needed  for  initial  conditions  for  direct  numerical  simulation  of  transition. 
Grid  adaptation  based  on  refinement  yields  better  than  power  law  convergence  in 
the  mean  flow  error. 

The  structure  of  eigenfunctions  in  different  branches  of  the  global  spectrum  is  inves¬ 
tigated.  The  adjoint  problem  is  solved  simultaneously  at  negligible  computational 
cost  and  the  structure  of  the  adjoint  eigenfunctions  analyzed.  It  was  observed  that 
the  normalized  adjoint  pressure  eigenfunction  is  typically  four  orders  of  magnitude 
larger  than  the  pressure  eigenfunction,  indicating  that  the  flow  is  highly  sensitive 
to  mass  sources. 


f  Graduate  Student,  Dept,  of  Aeronautics  &  Astronautics 
J  Professor,  Dept,  of  Mechanical  &  Civil  Engineering. 
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2.  Introduction 


The  numerical  schemes  in  use  for  solving  compressible  linear  stability  equations  can 
be  broadly  classified  into  boundary  value  methods  (BVM)  and  initial  value  meth¬ 
ods  (IVM).  IVMs  require  good  guesses  of  the  eigenvalue.  Although  they  require 
less  memory,  there  is  a  risk  of  missing  some  modes.  On  the  other  hand,  BVMs 
[Malik  &  Orszag  (1981),  Malik,  Chuang  Sz  Hussaini  (1982),  Malik  (1990)]  require 
computationally  intensive  global  searches  and  efficient  iterative  local  searches.  Spec¬ 
tral  methods  are  not  a  good  choice  due  to  the  movement  of  the  critical  layer  from 
0.2^  0.9^,  where  6  is  the  boundary  layer  thickness,  as  the  Mach  number  increases 

from  0  to  10.  Asymptotic  boundary  conditions,  which  are  analytically  derived,  are 
used  in  the  free-stream.  Spurious  unstable  modes  are  eliminated  as  in  Malik  (1982). 
Sutherland’s  law  of  viscosity  is  used  for  the  functional  dependence  of  viscosity  on 
temperature.  Keller’s  box  scheme  (second  order)  is  used  for  comparison. 

The  growth  or  decay  of  infinitesimal  perturbations  superposed  on  laminar  solu¬ 
tions  of  the  Navier-Stokes  equations  is  the  subject  of  linear  stability  theory.  In  this 
work,  the  basic  equations  governing  the  linear  stability  of  parallel-flow  compressible 
boundary-layers  are  derived  by  linearizing  the  Navier-Stokes  equations  about  the 
laminar  flow.  These  perturbations  are  assumed  to  be  of  the  form 

u'{x,y,t)  =  (1.1) 

For  temporal  stability  analysis,  a,  the  streamwise  wavenumber,  is  fixed  and  real  and 
O’,  the  frequency,  is  complex;  for  spatial  analysis,  a  is  complex  and  u  is  fixed  and 
real.  In  temporal  analysis,  u  =  Ur+iiOi,  Ur  is  the  frequency  and  Ui  is  the  growth  rate 
of  the  perturbation.  These  infinitesimal  disturbances  are  imposed  on  the  compress¬ 
ible  Navier-Stokes  equations  linearized  about  a  laminar  boundary  layer  solution.  If 
it  is  assumed  that  the  mean  flow  is  locally  parallel,  a  set  of  five  ordinary  differential 
equations  is  obtained.  Of  these,  three  are  the  second  order  momentum  equations, 
one  is  the  second  order  energy  equation,  and  one  is  the  first  order  continuity  equa¬ 
tion;  thus  the  complete  system  is  of  ninth  order.  For  reviews  of  boundary-layer 
stability  theorj^,  see  Reshotko  (1976),  Mack  (1984)  or  Nayfeh  (1989). 
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LINSTAB  is  a  compressible  linear  stability  analysis  code  for  two-dimensional  bound¬ 
ary  layers.  It  uses  an  iterative  finite-difference  method  to  compute  the  most  unstable 
eigenvalue  and  requires  an  accurate  estimate  of  the  most-unstable  eigenvalue.  A  lo 
cal  eigenvalue  search  procedure  improves  the  accuracy  of  the  eigenvalue  and  yields 
the  eigenfunctions.  A  global  procedure  was  developed;  it  may  be  used  when  no 
estimate  of  the  most  unstable  eigenvalue  is  available. 

The  solution  of  the  adjoint  problem  generally  identifies  regions  which  are  important 
for  the  placement  of  the  forcing.  The  adjoint  eigensolution  defines  the  efficiency  with 
which  a  particular  forcing  excites  the  eigensolution. 

3.  Mean  Flow 


3.1  Governing  equations 

The  Navier- Stokes  equations  governing  the  flow  of  a  viscous  compressible  ideal  gas 
are 


(3.1) 

(3.2) 

(3.3) 

(3.4) 

where  q  is  the  velocity  vector,  p  the  density,  p  the  pressure,  9  the  temperature,  R 
the  Reynolds  number,  Ec  (-  ul/{CpTe))  the  Eckert  number,  a  the  Prandtl  number, 
p,  the  coefficient  of  viscosity,  and  ^  the  viscous  dissipation  given  by 


l  +  U.v),  =v. 
|  +  v.(„)  =  o, 


P 


dR 


R' 


dt 


+  {q.W)d 


1 

P  =  pO 


^V.{pVe)  d-  Ec  <1  ^  +  (g.V)p  + 
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The  quantities  have  been  nondimensionalized  with  respect  to  the  frce-strcam  values. 
For  simplicity,  the  Stokes  approximation  of  zero  bulk  viscosity  has  been  assumed 
(A/,,  =  -2/3). 


The  mean  flow  is  taken  to  be  a  similarity  solution  of  the  boundary  layer  on  a  flat 
plate  with  no  pressure  gradient,  obtained  from  the  boundary  layer  approximation 
of  the  Navier-Stokes  equations  (3.1  -  3.5).  The  equations  from  which  this  solution 
can  be  derived  are  obtained  with  the  aid  of  the  Mangler-Levy-Lees  transformation 
(White  1991,  p.  534)  for  two-dimensional  flow 

=  p^HeUedx, 
dr]  =  pWe/(20’^^  dy- 

If  constant  properties  are  assumed,  Eqs.  (3.1  -  3.4)  reduce  to 

(c/" )'  +  ff"  =  0  (3.8) 

(aig'  +  a2f'f")'  +  fg'  =  0  (3.9) 

where  a  prime  denotes  differentiation  with  respect  to  rj. 

The  boundary  conditions  become 


(3.6) 

(3.7) 


m  -  /u.,  /'(O)  =  0,  lim  =  1,  (3.10) 

g{^)-gu>,  lim5r(77)  =  l.  (3.11) 

rj-^oo 

In  these  equations: 


/'  =  u/ue, 

g  =  H/He, 


c  =  pplpefi^, 
a\  =  c/<7, 


c 


and  H  is  the  enthalpy,  7  the  ratio  of  specific  heats,  p,.^  pe-,He  the  edge  values  of 
the  velocity,  density,  viscosity  and  enthalpy,  and  M  the  edge  Mach  number  defined 


as  M 


— The  Prandtl  number  a  is  defined  as  cr  =  where  Cp  is  the 
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specific  heat  at  constant  pressure  and  is  assumed  to  be  constant.  The  viscosity  is 
assumed  to  be  given  by  the  Sutherland  formula 


=  1.46  X  10"® 


jrl/2 

1  +  110.3/T 


N.sm 


The  thermal  conductivity  k  is  computed  using  a  —  0.7.  Atlternatively,  the  profiles 
could  be  provided  numerically  by  the  user.  For  example,  they  could  be  obtained 
from  a  two-dimensional  boundary  layer  calculation. 


Compressible  Linear  Stability  Equations 

Now  we  describe  the  procedures  for  temporal  stability  analysis.  The  procedures  for 
spatial  analysis  are  the  same  but  the  roles  of  u)  and  a.  are  exchanged.  Following  Eqs. 
(3.1)-(3.5),  let  Q,  P,  T  and  pm  be  the  velocity,  pressure,  temperature  and  density  of 
the  steady  mean  flow,  and  q',p\T',  and  p'  the  velocity,  pressure,  temperature  and 
density  of  the  disturbance.  Then  q  =  Q  +  q',P  =  P+P^  ^  =  T  +  T',  p  =  Pm  +  p'  and 
Q,  P,  T,  Pm  satisfy  the  Eqs.  (3.1)-(3.5).  Subtracting  the  equations  for  the  mean 
flow  from  the  full  equations,  linearizing  around  the  mean  flow,  assuming  the  mean 
flow  to  be  locally  parallel  (Q  =  {U{y),0,W{y)),  no  pressure  gradient  across  the 
boundary  layer  (p  =  1  and  pm  =  l/T),  and  the  disturbance  to  have  the  form 

«(y)  \ 

v(y)  j  exp  [i(Q;a;  +  /3z  —  ut)] , 
w\y)) 

p'  =  p(y)  exp  [i(Q;x  +  j3z  —  u>t)] 

T'  =  T(y)exp  [i(aa:  +  ^z  —  ivt)] .  (3-12) 


to  obtain  the  following  equations. 


+  PW  -  .)  +  +  /?^)]  +  (if  -  ^|a)  f,  +  iop 

h  1  f  5  /  ,dU\  ~  dU  dU 

+  ^po/?«'  ;j  I  aj,  (/<  a J  r  + 

dh  ,dU  df  a^-'j 

~  ^)  +  ■*■  ^  ■*■  ^ 

f  du  ^dw\  dn  dv  .  d^v\ 

+  1^/^  (  0-- — h  +  h^^—^  r  =  O’ 


(3.13) 


5j/  '  "  ay 


dy 


(3.14) 


ii(acr  +  0W  -  u,)  +  |(q^  +  +  (^i^  -  ^|i/5)  f.  +  \DP 

h  '^\df,dW\~  dW  dw 

R\dy\^  dy  du 


.,  ^S5  ,airaf  a^ii'l  „ 

1  ^  1. 
i-(aw  +  ^iB)  +  i(Q^7  +  -a;)p+  =  0’ 


(3.15) 


(3.16) 


\^{aU  +  m-u) 


(7-i)m2  ,)/aa' 


r  1  ar  21(7-1)11/2  /  ou  ^aii^M 

-  i(7  -  l)M\aU  +  -  u;)p 

2(7  -  i)m2  j^T  ayar  a^r 

^  ay  ay  ay  ay  y  Ra  dy  dy  ^  dy'^ 


0.  (3.17) 
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where  7  is  the  ratio  of  the  specific  heats,  cr  the  Prandtl  number,  R  theHeynolds 
number  based  on  the  displacement  thickness  6%  M  the  Mach  number,  and  Ij  = 

j  +  A/ IX. 

With  the  above  asumptions,  the  equation  of  state  (3.4)  simplifies  to 

p  =  (3.18) 

which  was  used  to  eliminate  density  p  from  (3.13)-(3.17). 

The  boundary  conditions  are 

u  =  v  =  w=^T  =  0  at  y  =  0,  (3.19) 

u,v,w,T  ^  0  as  y CO.  (3.20) 

They  are  the  linear  stability  equations  for  the  compressible  parallel  flow. 

3.2  Numerical  solution  schemes 

To  solve  Eqs.  (3.8)-(3.9),  Keller’s  Box  Method  is  used.  This  method  is  described  in 
great  detail  by  Cebecci  and  Smith  (1974).  We  will  apply  it  to  the  two-dimensional 
flow. 

We  first  write  Eqs.  (3.8)  and  (3.9)  as  a  first-order  system  of  ordinary  differential 
equations.  For  that  purpose,  we  introduce  new  dependent  variables  u(r))  and  v{rj) 
so  that  the  momentum  equation  (3.8)  can  be  written 

f  =u  (3.21a) 

u  =  V  (3.216) 

{cv)'  +  fv  =  0  (3.21c) 

In  terms  of  these  variables,  the  boundary  conditions  become 

/(0)  =  /^,  u(0)  =  0,  ^l^u(ry)  =  l,  (3.22) 
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and,  introducing  a  new  function  G{r]),  the  energy  equation  (3.9)  can  be  written 


(oiG)'  +  fG  +  {a2vvy  =  0 
g'  =  G 


along  with  boundary  conditions 


(3.23a) 

{3.23b) 


ff(0)  =  Qw  or  G(0)  =  Gu.,  lim  g{ri)  =  1,  (3.24) 

t;— »oo 

3.2.1  Second-order  difference  scheme 

The  discretization  is  based  on  Keller’s  Box  Method. 

Momentum  equation 

We  denote  {f,u,v)  at  (gj)  by  {fj,Uj.,Vj).  The  discretized  form  of  Eqs.  (3.21)  is 
obtained  by  integrating  each  equation  from  to  gj  and  using  the  trapezoid  rule 
to  approximate  any  integral: 


fj  fj-i  +  '^j-i)  —  0) 

(3.25a) 

icv)j  -  (cy)j-i  +  {{fv)j  +  (/v)_,_i)  =  0, 

(3.255) 

Uj+i  -  Vj  -  +Vj)  =  0, 

(3.25c) 

The  boundary  conditions  are 

/!=/«.,  wi=0,  uj  =  l,  (3.26) 

and  hj  —  gj  —  gj-i. 

A  Newton  iteration  method  is  used  to  solve  the  non-linear  system  (3.25)  along  with 
boundary  conditions  (3.26).  If  the  superscript  i  denotes  the  value  of  the  variable  at 
the  z'th  Newton  iteration,  then 

/r” = f'3 + 

^  uf + suf 
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Linearizing  the  system  (3.25)  about  the  solution  at  iteration  i  by  dropping  the  terms 
quadratic  in  ^-quantities,  we  obtain  the  linear  system 


^fj  “■  ^/j-i  “  +  <5wj-i)  —  rj, 

CjSvj  Cj—iSvj—i  Sj, 

6uj+i  —  Suj  —  ^hj+i{Svj^i  +  6vj)  =  tj, 

where 

rj  =  -  (/i  -  fj-i  - 

•Sj  =  -  -  Cj-iVj-i  +  \hj  {fjVj  +  fj-iVj-i)) 

tj  =  -  (uj+i  -Uj-  lhj+i(vj+i  -t-  Vj)) 

The  system  (3.27)  is  linear  and  has  block  tridiagonal  structure  and 
written  as 

+  AjSj  -f  Bj6j_i  -  Tj, 


where 


■  ■ 

Suj 

’  -i 

8vj 

- 1 

_ 1 

and,  for  j  =  2, 3, ...,  J  —  1, 


'0 

0 

0 

0 

0 

0 

[0 

1 

?+l . 

’1 

— 

0 

Aj  = 

0 

0 

C; 

h, 

- 

.0 

-1 

-\h 

j+l . 

■-1 

-^hj 

0 

0 

0 

Cj-1 

hi 

0 

0 

0  ^ 

The  boundary  conditions  are 

Sfi  =  0,  6ui  =  0  and  8uj  =  0 


(3.27a) 

(3.276) 

(3.27c) 


may  be 
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so  that,  at  the  wall, 


’l 

0 

0 

'0 

0 

0 

■  0  ' 

- 

0 

1 

0 

,  C,  = 

0 

0 

0 

>  Ill  = 

0 

0 

-1 

-1^2. 

0 

1 

- 1 

(N 

1 

^1  _ 

and,  at  the  free-stream. 


-\hj 

0  ■ 

■-1 

-\hj 

0 

rj 

Aj  = 

0 

0 

CJ 

hj 

to 

II 

0 

0 

cj-l 

hj 

>  r.j  = 

0 

1 

0 

0 

0 

0 

0  _ 

This  system  is  solved  by  means  of  block-tridiagonal  factorization  as  described  in 
Cebeci  and  Smith  (1974,  §7).  The  quantities  {fj,Vj,Vj)  are  updated  with 
(Sfj,6uj,6vj)  and  the  process  is  repeated  until  convergence  within  preassigned  tol¬ 
erance  (10“^°)  is  achieved. 


Energy  equation 


The  energy  equation  is  coupled  to  the  momentum  equation  through  the  fluid  prop¬ 
erties  like  viscosity  and  thermal  conductivity.  However,  once  these  parameters  and 
f,u  and  V  are  known,  the  energy  equation  becomes  linear. 

We  assume  initial  enthalpy  and  velocity  profiles  and  compute  the  fluid  properties 
and  then  solve  momentum  equations  followed  by  the  energy  equation  iterating  until 
the  convergence  criterion  is  met. 


Eqs.  (3.23)  are  linear  and  their  discretized  form  is 
(oiG)j  — 


hj 


+  +  fj-iGj-i)  -  tj, 

9j+i  -  9j  +  Gj)  =  0, 


(3.280) 

(3.281) 


along  with  boundary  conditions 

Gi  =  Gw  =  0  and 
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9J  =  1. 


As  for  the  momentum  equations,  Eqs.  (3.28)  are  written  in  a  block-tridiagonal  form 


where 


with 


4- 


93 

^3 

[GiJ 

’  ^3  - 

0 

tj  =  -—((a2uv)j  +  {a2uv)j-i) 


and,  for  j  =  2, 3, ...,  J  —  1, 


Ci 


■^3  ~ 


0  0 

1  J 

0  + 

-1  J 

n  _ (“ilj-i  1  i  f . 

hi 


0  0 

The  boundary  conditions  give  rise  to  special  cases:  At  the  wall, 

0  0 

2^2 


0  1 

-1  -\h2 


,  Ci  = 


1  -Ihi 


,  ti  = 


wile,  at  the  free  stream, 

(°i)j  I  1 


Aj  = 


0  ^  + 
1  0 


,  Bj  = 


(ai)  J-i  j_  1 


hj  '  2 

0  0 


+  2fj-^ 


5  LJ 


r  T  = 


tj 

1 


This  system  is  solved  by  means  of  block-tridiagonal  factorization. 
3.2.2  Fourth- order  compact  difference  scheme 


A  fourth-order  accurate  two-point  scheme  can  be  derived  by  means  of  the  Euler- 
Maclaurin  formula 


_  ,[,‘-1  =  h  (^  +  +  0(hl),  (3,29) 

2  \  dy  dy  J  12  \  dy^ 


dy"^ 
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where  and  =  yu  -  Vk-i- 


This  high-order  numerical  method  was  first  used  for  boundary  layers  by  Wornom 
(1977)  who  demonstrated  its  efficiency  compared  to  other  fourth-order  methods, 
especially  when  non-regular  meshes  are  used.  Iyer  and  Harris  (1989,  1990)  present 
computations  using  this  compact  scheme  for  three-dimensional  compressible  bound¬ 
ary  layers. 

The  strong  point  of  this  scheme  is  that  only  two  data  points  are  required  and  grid 
smoothness  is  not  crucial  for  high  accuracy.  Second  order  local  computations  are 
sensitive  to  stretching  functions. 


Momentv,m  equation 


We  apply  formula  (3.29)  to  Eejs.  (3.21)  with 


h 

Vk 


Uk 


so  that 


where 


— —  =  and 

dy 


d^xijk 

dy"^ 


Ak  - 

0 

_  £* 

Ck 

0 

_C_[ 

Ck 

r 

0 

,  = 

0 

fkVk  , 

+2 

Ck  J 

1 - 

1 

0 

1 

0 

—  2Jl 

_£k 

0  . 

Ck 

Ck 

0 

0 

r 

o 

o 

o 

II 

_  ILL 

Ck 

_£*. 

Ck 

0 

II 

1 

0  0  0 

? 

0 

0 

0 

0  1  0 

0 

1 

0  ■ 

0 

0 

O' 

to 

II 

fkVk  1 

^//  /  /  \  ^ 
+2 

Ck  '  \‘^l'  J 

—  ILL 
Ck 

0 

0 

c', 

0 

1  0 

0 

0  . 

_ ^ 

Ck 

_ LL 

Ck 

0 

5 
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so  that 


fk  fk  —  1 
Vk  -  Vk-1 
Uk —  Uk+1 

_h±i  *‘)  +  %  .  (3.30) 

2 

These  equations  can  be  written  in  the  block-tridiagonal  form 

+  =  0  (3.31) 


=  h  ^  (B?®‘  -  B+.J* 


fe-1' 


where 


P,  =  I-'}±At+  'if- a:  +  +  %.B,-, 

«i  =  -  (/+  +  ^4-1  +  ^Bti)  . 

Rk  =  -(r-  fi-Ai 


+ 


/,2  ' 
^+1  R- 


(3.32a) 

(3.326) 

(3.32c) 


with 


J+  = 


‘1 

0 

o' 

'0 

0 

0' 

0 

1 

0 

and 

r  = 

0 

0 

0 

0 

0 

0 

0 

0 

1 

In  order  to  solve  this  system  by  Newton  iteration,  we  need  to  linearize  Eq.  (3.31). 
If  i  denotes  the  number  of  the  Newton  iteration,  then 


where 


Ak  = 


Bk  = 


Vk 

Ck 

0 


0 

ifk  +  cj) 

Ck 


1 

0 

0 


0  ^  2  ° 

^(34+2/0  if(34  +  A)  - +  2  (I)  t 

(/fc+c'fe)  n 


—  Uk. 
Ck 


Ck 
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Again  we  separate  Ak  and 


At 
Bt  = 
b;  = 


0 

0 

r 

'0 

0 

o' 

(fk+c'k) 

0 

0 

0 

0 

Ck 

Ck 

0 

0 

0 

0 

1 

0 

0  1 
7f(3c',  +  2/,)  +  + 

0  *  0 


0  0  O' 

0  0,0 

Vk  (fk  +  c't)  Q 

Ck  Ck 


0  1 


ILL 

Ck 

0 


so  that  we  now  have  the  linear  equation 


+  PkS<i>^  +  =  Hk  (3.33) 

where  Pk^  Qk  and  Rk  are  defined  from  Ak  and  Bk  as  in  (3.32)  and 

Hk  =  -  +  Pk^^  +  Qk^^-^) . 

The  boundary  conditions  (3.22)  require  the  following  modifications. 


At  the  wall, 


'1 

0 

O' 

'0 

0 

0' 

0 

Pi  = 

0 

0 

1 

,  — 

0 

0 

0 

,  H,  = 

0 

- 

same  as  A-  /  1 

- 

. 

same  as  A-  1 

same 

while,  at  the  free  stream, 

same  as  k  ^  J 

1 

- 

same  ask  ^  J 

same 

Pj^ 

same  as  k  ^  J 

,  Qj  = 

same  as  k  ^  J 

,  Hj  = 

same 

0 

0 

1 

0 

0 

0_ 

0 

This  system  is  solved  by  means  of  a  block- tri diagonal  factorization.  is  updated 
with  8'^^  and  the  process  is  repeated  until  convergence  within  preassigned  tolerance 
is  achieved. 
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Energy  equation 


For  the  energy  equation  (3.23)  we  follow  the  same  steps  as  above  (the  linearization 
is  unnecessary  as  Eqs.  (3.23)  are  already  linear): 


As  Eqs. (3. 24)  are  not  homogeneous,  we  write 


I 

——-Ak'^^'  +  ^k  and 
dy 


d^ 

dy'^ 


=  Bk-^^  + 


d^k 

dy 
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As  for  the  momentum  equations,  the  system  can  be  written  in  the  block  -  tridiagonal 
form 


where 


P,  =  I-^At  + 


hk  .+  ,  hk+i 


AT  +  -^BT  + 


k  D+  I  ^  +  1  p- 


2  -  2  ■"12^^  '  12 


k  ’ 


0,.  =  -  (n  +  . 


Rk  = 


hk+i 

-^^k+i 


at,.  + 


12 


*:  +  l  j  > 


(3.34) 

(3.35a) 

(3.355) 

(3.35c) 


and 


‘  2  ^  ‘  *■'>'  12  ds  dy 


2  ^  ^  12  \  dy  dy 

The  boundary  conditions  (3.24)  require  the  following  special  cases.  At  the  wall, 


Pi  = 


same  as  A’  ^  1 

0 


1 


,  Ri  — 


while,  at  the  free  stream. 


Pj 


same  as  k  ^  J 


,  Q. 


same  as  A’  ^  1 

0  0 


0  0 
same  as  k  ^  J 


= 


,  Hj 


same 

0 


1 

same 


This  system  is  solved  by  means  of  a  block-tridiagonal  factorisation. 


4.  Global  Method 


When  no  guess  of  the  most-unstable  eigenvalue  is  available,  LINSTAB  uses  a  global 
method  that  computes  the  entire  eigenvalue  spectrum.  The  second-order  finite  dif¬ 
ferenced  compressible  stability  equations  (3.13)-(3.17)  can  be  reformulated  as  a 
matrix  eigenvalue  problem 

A$=a;B$,  (4.1) 
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where  u  is  the  eigenvalue  and  $  the  discrete  representation  of  the  eigenfunctions 


'au  +  I3w' 

V 


L  aw  —  ^u] 


(4.2) 


The  eigenvalues  are  the  roots  of  the  determinant  equation 


Det\A-uB\  =  0.  (4.3) 

This  is  a  standard  matrix  eigenvalue  problem  and  is  solved  using  the  complex 
‘shifted’  LR  method  described  in  Wilkinson  (1965). 


The  most  unstable  eigenvalue  is  the  one  with  largest  imaginary  part.  The  global 
eigenvalue  search  is  formulated  so  that  no  growing  spurious  modes  (not  physically 
relevant)  are  generated,  see  Malik  (1982).  The  spurious  modes  are  eliminated  by 
using  a  finite-difference  scheme  for  the  eigenvalue  analysis  that  is  consistent  with 
the  scheme  for  the  initial-value  problem.  When  the  eigenvalue  problem  (4.3)  is 
solved  using  the  complex  LR  algorithm,  the  storage  requirements  are  0{K‘^)  while 
the  computational  work  is  0(/\^)  where  K  =  5N  for  the  eighth  order  system.  The 
global  method  is  computationally  expensive  and  should  be  used  only  when  no  guess 
of  the  eigenvalue  is  available. 

All  discrete  eigenvalues  obtained  using  complex  ‘shifted’  LR  method.  The  linear 
disturbance  equations  constitute  an  eigenvalue  problem  which  yields  the  complex 
dispersion  relation  w  =  uj(^a,l3')  and  the  construction  of  this  relation  is  the  major 
concern  of  linear  stability  analysis. 

5.  Local  Method 

When  a  guess  for  the  most-unstable  eigenvalue  is  available,  it  can  be  improved  by 
a  local  method  which  also  computes  the  corresponding  eigenfunction.  In  this  sec¬ 
tion,  we  present  two  ways  of  numerically  solving  the  compressible  stability  problem; 
inverse  Rayleigh  iteration  and  Newton  iteration. 
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5.1  Inverse  Rayleigh  iteration 


This  procedure  was  used  in  the  original  version  of  LINSTAB,  and  its  theory  is 
presented  in  Wilkinson  (1965).  Generalization  of  this  procedure  to  the  compressible 
stability  problem  (3.13-3.17)  results  in  the  following  algorithm 


(A  - 

(5.1) 

(A  - 

(5.2) 

(5.3) 

The  error  satisfies  cja  +  i  —  ^  —  tu)^).  The  iteration  cycle  is  started  with 

the  guessed  eigenvalue  produced  bj^  the  global  method,  loq,  and  an  assumed  but 
arbitrary  smooth  profile  for  the  eigenfunction  $(0)  and  its  adjoint  ^'(0).  The  algo¬ 
rithm  converges  cubically  for  the  eigenvalue,  but  the  eigenfunction  converges  at  the 
square  root  of  this  rate,  as  shown  bj"  Hackbusch  (1985). 

5.£  Newton  iteration 


Another  method  of  local  eigenvalue  search  to  improve  the  accuracy  of  the  most 
unstable  eigenvalue  is  Newton  iteration,  which  yields  the  same  level  of  accuracy 
for  the  eigenfunction  and  the  eigenvalue.  In  symbolic  form,  the  system  of  ordinary 
differential  equations  satisfied  by  the  linear  disturbances  can  be  written 


=  H, 


(5.4) 


where  $  is  defined  by  eqn.  (4.2)  and  H  =  0.  The  boundary  conditions  for  Eq.  (5.4) 
are 


j/  =  0;  (f>;i  =  (1)2  =  =  (f>5  =  0 

y  ^  oo; 


(5.5) 


For  the  local  eigenvalue  problem,  Eqs.  (5.4)  are  a  block-tridiagonal  system  which 
is  solved  using  LU  factorization.  As  Eq.  (5.4)  is  homogeneous,  in  order  to  avoid 
a  trivial  solution,  one  inhomogeneous  boundary  conditions  is  imposed  at  the  wall. 
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Specifically,  as  proposed  by  Malik  (1990),  the  boundary  condition  <^i(0)'^  0  is  re¬ 
placed  by  <j)3{0)  =  1.  This  is  equivalent  to  normalizing  the  eigenfunction  so  that 
the  pressure  perturbation  at  the  wall  is  unity.  Since  the  pressure  does  not  vanish 
at  the  wall,  this  condition  is  appropriate;  see  Malik  (1990)  for  other  possible  nor¬ 
malizations.  A  non-trivial  solution  may  now  be  obtained  if  u;  =  wq,  the  correct 
eigenvalue.  Newton’s  method  is  used  to  iterate  on  w  until  the  missing  boundary 
condition  =  0  is  satisfied.  After  a  solution  $  is  obtained  using  the  estimated 

eigenvalue  coq  ,  the  correction  Ao;  is  determined  from  the  linearized  equation 

^,(0)+^/lu,  =  0.  (5.6) 


where  ^i(O) 
solving 


is  known  from  the  solution  $  just  computed;  d(j)i{0)l du)  is  obtained  by 


9$  dL 


(5.7) 


The  process  is  repeated  until  (f>\(0)  vanishes  within  a  preassigned  tolerance. 


Boundary  conditions 


Homogeneous  boundary  conditions  are  used  except  for  the  pressure.  The  tempera¬ 
ture  perturbation  is  set  to  zero  at  the  solid  boundary.  This  is  a  reasonable  assump¬ 
tion  since  high  frequency  disturbances  do  not  penetrate  into  the  wall  due  to  thermal 
inertia.  In  other  words,  the  wall  appears  insulated  on  the  time  scale  of  mean  flow 
but  not  on  the  short  time  scales  of  the  disturbances. 


5.2.1  Second-order  difference  scheme 

The  second-order  scheme  used  to  discretize  the  linear  disturbances  equations  (5.4) 
together  with  boundary  conditions  (5.5)  is  based  on  the  following  system  of  ordinary 
differential  equations 

(AI>2  +  HD -b  C)$  =  0.  (5.8) 

Here  D  =  d/dy,  while  A  is  given  as 
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n 


)  ij 

and  B  and  C  are  5x5  matrices  whose  nonzero  elements  are  given  in  Appendix  I. 


The  system  is  discretized  using  second-order  central  differencing  and  leads  to  a 
block-tridiagonal  system  which  is  solved  using  a  block  LU  elimination  (Thomas 
algorithm). 

5.2.2  Fourth- order  compact  difference  scheme 

To  reach  fourth-order  accuracy  we  use  the  two-point  compact  difference  scheme 
presented  in  Section  2.2.2., 


In 

2  V  dy  ^  dy  ) 


12  V  dy2  ^y2  J 


+  o(/4), 


(5.9) 


where  =  ^'(yA.)  and  h/,  =  y^  -  yr-i. 


Using  the  continuity  equation,  the  second-order  normal  momentum  equation  may 
be  reduced  to  a  first-order  equation  for  pressure.  Thus  the  linear  stability  equations 
(5.8)  may  be  rewritten  as  a  system  of  eight  first-order  equations 


dy 


i  =  1,2,  ...,8, 


where 

=  au  -I-  I3w,  ip2  =  ip3  =  V, 

dy 

=T,  ipe  =  ipr  =  aw  -  ^u, 

dy 

with  corresponding  boundary  conditions 


i’i  =  p, 


(5.10) 


y  0;  V’l  =  V’S  =  V’5  =  •07  =  0 
y  ^  oo;  ip^^tp3,ip5,tp-!  ^  0. 


(5.11) 
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The  coefficients  aij  are  given  in  Appendix  II.  In  order  to  apply  the  compact  scheme 
to  Eq.  (6.4),  we  also  need 


d?xl) 


dy‘ 


i=i 


which  is  obtained  by  differentiating  Eq.  (4.10).  Here 

8 


,  daij 


+  ^  O-ilO-lj 


1-1 


(5.12) 


Thus  Eq.  (5.9)  becomes 


hk 


2  ^  1 


hi 


8 


k 

22  /  ^  ^  j 


8 


hi 


j-i 


12 


j=i 


=  0  (5.13) 


or 

=  0.  (5.14) 

Eqs.  (5.14)  along  with  boundary  conditions  (5.11)  are  then  written  in  the  block- 
tridiagonal  form 

+  Ak'^'' +  Ck^^'^^  =  H,  (5.15) 

where  Ak,  Bk,  Ck  are  8x8  matrices  defined  below  and  i?  is  a  8  x  1  null  matrix. 


Ak 


Bk 


Ck 


last  4  rows  of  Rk 

first  4  rows  of  Sk+i 
last  4  rows  of  Sk 

0 

0 

first  4  rows  of  Rk-f  i 


A:  =  2,iV-l, 

k  =  2,iV, 

fc  =  l,iV-l, 


(5.16) 


(5.17) 


(5.18) 
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>ll  - 


E 


first  4  rows  of  S2 


An 


last  4  rows  of  Rk 


F 


(5.19) 


(5.20) 


where  E  and  F  are  4x8  matrices  representing  the  bottom  and  top  boundary 
conditions  (4.11).  Thus 


£^'(0)  =  0  and  F^'(oo)  =  0, 


(5.21) 


where 

(5.22) 

and  the  elements  of  F  are  determined  from  the  asymptotic  behavior  of  as  y  — >  00. 
For  the  details  of  the  computation  of  F,  see  Malik  (1982  -  App.  II)  and  Mack  (19C5, 
pp  266-271).  As  discussed  above,  the  nonhomogeneous  boundary  conditions  needed 
to  avoid  a  trivial  solution  when  solving  system  (5.15)  is  obtained  by  setting 


0  0  1  0  0  0  0  0 

0  0  0  0  1  0  0  0 

.0  0  0  0  0  0  1  0. 


-0  0  0  1  0  0  0  0- 

0  0  1  0  0  0  0  0 

0  0  0  0  1  0  0  0’ 

.0  0  0  0  0  0  1  0. 


(5.23) 


and  H  =  (1,0,0,...,  0)^. 


The  block-tridiagonal  system  (5.15)  can  now  be  solved  by  block-LU  elimination. 


5.2.3  Acciiracy: 


For  a  second  order  method. 


<i^h  —  <f>ex  +  + 


4>2h  —  <I>€X  +  4:h^1pl  -1- 


23 


(j)^^  _  (‘i4>h  ^<i>2h)  truncation  error  ratio,  e  = 
For  a  fourth  order  method, 


(t>h  —  <t>ex  +  + . 

<l>2h  —  4>ex  +  + . 

(t>ex  =  and  truncation  error  ratio,  e  =  e x -t O ' '  ~ 


5.3.2  Group  Velocity  Computation 


If  (a,  /3)  is  the  wave-vector,  then  group  velocity  is  given  by 


du  du! 


The  compressible  stability  equations  for  three-dimensional  disturbances  are  in  the 
form 

L{a,l3,u;{a,/3))f  =0. 

After  the  above  equation  is  differentiated  with  respect  to  w,  taking  the  inner  product 
with  the  adjoint  xf  of  the  eigenfunction  $  gives  an  expression  for  the  group  velocity. 
Since  =  =  0 


du  _  (xf,  |^(f) 

This  method  of  computing  the  group  velocity  [Malik  1982]  is  more  efficient  than 
finite-differencing  viz.  (|^)j  =  velocity  is  obtained  at  less 

than  10%  of  the  cost  of  the  eigenvalue  search. 

6.  Grid  Adaptation 

The  grids  used  in  linear  stability  analysis  have  a  significant  effect  on  the  results. 
When  the  second-order  accurate  numerical  scheme  is  used,  the  grids  need  to  be  very 
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smooth.  In  the  original  version  of  LINSTAB,  exponentially  stretched  meshes  were 
used,  which  are  not  well  suited  for  grid  adaptation. 

By  considering  the  truncation  error  inherent  in  finite-difference  approximations, 
Vinokur  (1983)  proposed  grid  stretching  functions  based  on  the  inverse  hyperbolic 
sine.  Besides  being  smooth,  ‘Vinokur’  grids  allow  considerable  control  of  the  distri¬ 
bution  of  the  points,  which  is  not  the  case  with  exponentially  stretched  grids.  The 
slopes  at  the  wall  and  free-stream  can  be  specified.  It  is  also  possible  to  join  different 
Vinokur  grids  while  retaining  smoothness  at  the  junctions.  Finally,  the  number  of 
grid  points  does  not  affect  the  shape  of  the  grid.  Several  Vinokur  grids  were  tried 
in  the  local  computation.  The  resulting  change  in  the  eigenvalue  uj  is  an  indicator 
of  how  sensitive  these  computations  are  to  the  grid.  This  led  us  to  implement  an 
adaptive  grid  algorithm  controlled  by  the  error  in  the  eigenfunctions. 

Numerous  adaptive  grid  methods  are  available.  Most  can  be  divided  into  two  cate¬ 
gories:  displacement  methods  and  refinement  methods.  The  first  type  used  a  fixed 
number  of  mesh  points  and  the  adaptation  consists  of  moving  the  mesh  points 
from  low-gradient  regions  to  high-gradient  regions.  The  second  type  starts  with 
a  coarse  mesh  and  adds  points  in  high-gradient  regions.  In  an  attempt  to  mini¬ 
mize  the  number  of  mesh-points,  we  first  tried  the  displacement  method  based  on 
a  error  equidistribution  variational  process  proposed  by  Eiseman  (1987).  Althoiigh 
promising  at  low  Mach  number,  this  method  was  not  able  to  handle  the  very  steep 
gradients  in  the  hypersonic  second  mode  eigenfunctions,  especially  the  temperature 
eigenfunction.  The  main  cause  was  the  inability  of  the  algorithm  to  maintain  the 
proper  smoothness  of  the  grid. 

It  was  then  decided  to  develop  a  grid-refinement  method  that  would  maintain  the 
required  mesh  smoothness.  The  algorithm  is  defined  by  the  following  steps  : 

(i)  The  initial  grids,  Go  and  Gi,  are  Vinokur  grids  of  41  and  81  points; 

(ii)  The  eigenfunctions  computed  on  G,-  and  G,_i  are  compared  and  estimates  of 
the  truncation  error  are  constructed.  Refinement  intervals  are  introduced  where  the 
truncation  error  is  greater  than  a  specified  e; 
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(in )  The  number  of  points  on  each  refinement  interval  is  doubled; 

(iv )  Smooth  connection  between  new  grids  and  the  old  ones  is  assured  by  a  data- 
passing  scheme; 

(v)  Repeat  steps  (ii )  -  (iv )  until  no  more  refinement  intervals  are  found. 

The  main  difficulty  was  the  choice  of  interpolation  method.  To  avoid  the  wiggles 
that  appear  with  B-spline  interpolation,  we  used  the  interpolation  method  based  on 
piecewise  cubic  polynomials  proposed  by  Akima  (1970).  The  slope  is  determined  us¬ 
ing  a  second-order  geometric  rule  which  leads  to  very  “natural  looking’  and  smooth 
grids.  The  junctions  between  the  old  parts  of  the  grid  and  the  new  parts  generated 
with  Akima’s  method  were  made  from  fifth-order  B-splines. 

When  using  a  fourth-order  accurate  method,  the  smoothness  requirement  is  less 
severe  than  for  second-order  methods,  and  thus  Akima  interpolation  can  be  used 
everywhere. 

This  new  method  was  first  tested  on  a  Mach  0.5  case  for  which  the  eigenfunctions 
are  smooth.  The  iteration  process  converged  rapidly  (in  less  than  10  iterations)  for 
e  «  10“®.  For  lower  values  of  e,  the  first  10  iterations  simply  doubled  the  points 
everywhere  thus  producing  a  huge  number  of  grid  points.  This  demonstrated  the 
need  for  early  capture  of  the  significant  gradients.  Consequently,  to  obtain  very  low 
error  (say  e  «  10"'’'),  we  first  need  to  converge  for  an  intermediate  value  of  e  and 
then  restart  the  process  on  the  resultant  grid  with  a  smaller  value  of  e. 

Another  option  of  the  grid  adaptation  routine  permits  reduction  of  the  number  of 
points  when  a  new  structure  has  been  captured,  before  continuing  the  adaptation 
process  with  a  lower  value  of  e.  This  allows  one  to  keep  the  number  of  grid  points 
as  low  as  possible  while  building  some  structure  in  the  grid.  , 

7.  Results  Conclusions 

Fourth— order  accuracy  is  obtained  both  for  the  meanflow  and  the  most  unstable 
eigenfunction  for  the  temporal  stability  problem.  The  advantage  of  this  scheme  is 
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the  compact  stencil  which  requires  only  two  data  points  for  fourth-order  accuracy. 
It  is  also  quite  efficient  compared  to  second-order  methods  requiring  fewer  grid 
points  to  resolve  the  flow  accurately.  These  computations  are  not  sensitive  to  the 
type  of  stretching  functions  like  second-order  computations. 

Analysis  of  adjoint  pressure  eigenfunction  indicates  that  the  flow  is  sensitive  to  mass 
sources.  The  structure  of  the  eigenfunctions  is  qualitatively  different  in  different 
branches  of  the  global  eigenvalue  spectrum.  However,  the  modes  investigated  have 
a  localized  structure  close  to  the  wall. 

Locally  refining  the  grid  using  error  estimates  is  a  useful  tool  for  ensuring  accuracy 
of  the  eigenfunctions.  Grid  adaptation  based  on  refinement  shows  better  than  power 
law  convergence  in  the  error. 

The  mean  flow  profiles  for  the  compressible  flat  plate  with  an  adiabatic  wall  were 
compared  to  the  ones  of  Van  Driest  (1952).  The  quantities  compared  are  the  dimen¬ 
sionless  boundary-layer  thickness  and  the  mean  skin-friction  coefficient 

Cfy/Rcx-  Fig.  1  shows  the  boundary  layer  thickness  quadratic  dependence  on  Mach 
number  for  adiabatic  wall  conditions  with  fixed  Reynolds  and  Prandtl  numbers  and 
edge  temperature.  The  computed  boundary  layer  thickness  agrees  quite  well  with 
the  Van  Driest  curve. 

Fig.  2  shows  the  effect  of  Mach  number  on  non-dimensional  skin-friction  coefficient. 
The  effect  of  increasing  Mach  number  is  to  decrease  the  skin  friction  coefficient  by 
about  a  factor  of  two  as  the  Mach  number  increases  from  0  to  20.  The  wall  shear 
stress  [tu-  =  increases. 

Fig.  3  shows  the  meanflowfor  M  =  4.5,  Re  —  8000,  Pr  =  0.7.  The  thermal  boundary 
layer  thickness  is  greater  than  velocity  boundary  layer  thickness  for  Prandtl  number 
less  than  unity.  The  adiabatic  wall  causes  a  bulge  in  the  temperature  profile  since 
there  is  no  heat  transfer  through  the  wall  and  viscous  dissipation  heats  the  flow  in 
the  near- wall  region. 

Fig.  4  shows  the  mean  enthalpy  profiles  for  M  =  4.5,  i?e  =  8000,  Pr  =  0.7  with 
isothermal  and  adiabatic  wall  conditions.  For  adiabatic  wall  conditions,  the  profile 
has  zero  slope  at  the  wall  and  a  bulge  in  the  boundary  layer  as  expected. 
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Fig.  5  shows  the  meanflow  convergence  history  for  M  =  4.5,  Re  =  8000.  The  relative 
error  is  e  =  maxj[|(l  -  where  Uex  is  the  exact  solution  of  the  discretized 

equations  at  grid  point  j  and  e  is  the  relative  error.  Typically  the  relative  error 
drops  one  order  of  magnitude  in  per  Newton  iteration. 

Fig.  6  shows  the  variation  of  mean  flow  error  with  number  of  grid  points  on  a 
logarithmic  scale  for  second  and  fourth  order  methods.  The  slopes  are  as  expected. 

Fig.  7  compares  the  efficiency  of  second  order  scheme  to  that  of  the  fourth  or¬ 
der  compact  scheme  Newton  iteration  (N-R)  for  the  mean  flow  calculation.  For 
Re=2000,  the  fourth-order  scheme  gives  reasonable  accuracy  with  N=41  while  for 
N  >  400  is  required  for  the  second  order  scheme. 

Fig.  8  shows  that  for  a  fourth  order  method,  the  ratio  of  errors  on  successive  grids 
(number  of  grid  points  =  250,  500,  1000)  higher  more  16,  indicating  fourth  order 

accuracy. 

Fig.  9  shows  the  error  distribution  of  mean  flow  enthalpy  ej  =  \{hj  —  hex)\  for 
M  =  4.5,  Re  =  8000.  The  error  decreases  by  a  factor  of  2^  for  every  grid  doubling 
showing  that  the  scheme  is  indeed  fourth  order  accurate.  The  error  curves  collapse 
when  scaled  appropriately.  The  dip  in  the  absolute  error  is  due  to  a  sign  change. 

Fig.  10  shows  that  the  global  eigenvalue  spectrum  has  two  branches.  The  number 
of  eigenvalues  increases  with  increasing  grid  resolution.  The  structure  of  eigenfunc¬ 
tions  is  different  in  each  branch.  The  eigenfunction  structure  in  the  two  branches  is 

investigated. 

Figs.  11  and  12  show  the  eigenvalue  structure  in  two  branches  of  the  global  eigen¬ 
value  spectrum.  In  the  mode  shown  in  fig.  11  for  A  =  (1.085,-0.555),  the  struc¬ 
ture  is  localized  close  to  the  wall  while  the  mode  in  figure  12  for  the  eigenvalue 
A  =  (0.585, 25.93)  is  highly  oscillatory. 

Fig.  13  shows  the  non-dimensionalized  pressure  and  its  adjoint  eigenfunction  struc¬ 
ture.  At  Re=8000  and  M=4.5,  the  pressure  adjoint  is  about  4  orders  of  magnitude 
higher  indicating  that  the  flow  is  sensitive  to  mass  sources  or  errors  equivalent  to 
mass  sources. 


28 


Fig.  14  shows  the  temperature  and  its  adjoint  eigenfunction  structure.  At 
Re=8000  and  M=4.5,  the  structure  is  quite  confined  close  to  the  wall. 

Fig.  15  gives  the  discretization  error  of  the  most  unstable  eigenvalue  of  the  sec¬ 
ond  and  fourth  order  schemes  using  local  search.  If  (A^,Aj)  are  the 

most  unstable  eigenvalues  of  the  discretized  equations  and  on  a  grid  jf,  then  £  = 

\/(A“-Aj)2  +  (A"-A;)2 

Fig.  16  shows  the  comparison  of  the  discretized  error  (accuracy)  in  the  most  unstable 
eigenvalue  computed  by  the  fourth  order  compact  Newton  iteration  for  exponential 
and  Vinokur  stretching  functions.  Exponential  stretching  appears  to  give  slightly 
higher  accuracy  for  the  same  number  of  grid  points. 

Fig.  17  shows  the  CPU  time  vs.  accuracy  for  inverse  Rayleigh  iteration  on  a  SPARC- 
10  station. 

Fig.  18  shows  the  CPU  time  for  global  and  local  calculations  on  C-90.  The  cost  of 
the  global  method  scales  as  while,  for  the  local  method,  it  increases  as  N,  where 
N  is  the  number  of  grid  points. 

Fig.  19  shows  that  the  real  part  of  group  velocity  increases  linearly  with  Mach 
number  while  the  imaginary  part  seems  to  decrease  with  M  for  Re  =  8000,  a  — 
2.25,/?  =  0. 

Fig.  20  shows  the  comparison  of  an  estimate  of  the  absolute  maximum  error  with 
and  without  refinement  for  the  mean  flow.  The  adaptive  grid  error  decreases  faster 
than  exponentially. 
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Appendix  I 

The  non-zero  coefficients  of  the  5  x  5  matrices  B  and  C  of  Eq.  (5.8)  are  given  below. 


Tj  1  dfl  rpl 

■Oil  - 

B\2  =  iC-^  — 

Bu  = 

B21  =  i(A  —  1)/ A 

■022  -  IldT-'- 
D  _  R 

■^23  - 

Bz2  =  1 

^41  =  2(7  -  l)M2(T(at7'  +  /3TT')/(«^  + 

D  _  2  dfl  rpl 

-O44  —  ^ 

Bi5  =  2(7  -  l)AfV(aW”  -  /?C')/(a^  +  /?") 

B54  =  i|f,(air'-/3i7') 

D  _  1  da  rpl 

•055  —  ^  ^-1 

Cn  =  -  [^(aU  +  /JTT  -  O))  +  A(a2  + 

C12  =  -  [^(oC'  +  0W')  +  +  /?") 

Cis=-f{a"+/3") 

C,4  =  (o^'  +/3W')S314T'  +(c.!7"  +/?»"')i5|. 

r<  _ ;  A— 2  1  dn  rpi 

^21  —  1  A  /i  dT-^ 

C22  =  -  [#s(aC^  +fiW-w)  +  (a"  +  /?2)/A 

C24  =  iift(oC'  +  /JW') 

C'ai  =  i 


54 


C'33  =i-/M^aU  +  [3W  -uj) 


C34  — 
C42  — 


-j.{aU  +  /3W-u) 

-  -  2i(7  -  l)APa(aU'  +  ^W) 


C43  =  ^(7  -  l)MHaU  +  /3W  -  u) 


C44 
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-cu) 
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+^-) 
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n  dT^-^  n  dT^ 

C52 

=  -ffiaW  -  fiV) 

C54 

=  -  -  pv) 

C55 

=  — 

+  /?TT  - 

■Uj) 

+  (0^ 

+  /3^) 

The  primed  quantities  are  the  derivatives  with  respect  to  the  boundary-layer  coor¬ 
dinate  y.  A  is  defined  as  2/3(/i2  +  2)  where  p2  is  the  ratio  of  the  second  coefficient 
of  viscosity  to  the  first. 


All  the  velocities  have  been  scaled  by  Ue,  the  streamwise  component  of  velocity  at 
the  edge  of  the  boundary  layer,  and  all  the  lengths  have  been  scaled  by  <5*,  the 
displacement  thickness.  The  resulting  Reynolds  and  Mach  numbers  are  then  given 

by 

R  = 

M  = 

where  pc,  Pe  and  Te  are  the  density,  viscosity  and  mean  temperature  in  the  free 
stream. 


These  coefficients  are  for  the  3D  flow.  In  the  code  W  is  set  to  zero. 
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Appendix  II 


The  non-zero  coefficients  of  Eq.  (5.10)  are  given  below. 
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«63  =  ^T'  -  2i(7  -  l)MV(aU'  +  /?W') 
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where  (  )  s  d/dy,  (  =  (aU  +  -  co),  x  =  ^3  +  X/y 
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